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ABSTRACT 


New  data  from  the  Seismological  Research  Observatory 
(SRO)  is  used  in  conjunction  with  a  non-linear  least 
squares  technique  to  invert  surface  wave  group  velocity 
data  for  the  shear  velocity  structure  of  the  South-China 
subplate.  A  multi-filter  analysis  was  used  to  produce 
group  velocities  for  the  fundamental  mode  Love  and 
Rayleigh  wave.  A  number  of  earthquakes  over  a  single  path 
are  used  in  order  to  deduce  a  measure  of  the  observational 
uncertainty  of  the  present  method.  Group  velocity  standard 
deviations  range  from  .03  km/sec  to  .15  km/sec  for  Rayleigh 
waves  and  .03  to  .1  km/sec  for  Love  waves  over  the  period 
range  of  10  to  65  seconds. 

Results  of  the  inversion  indicate  that  a  3  layer, 

40  km  thick  crust  is  an  adequate  model  consistent  with  the 
data  for  the  South-China  subplate.  Crustal  shear  veloci¬ 
ties  are  found  to  be  higher  and  upper  mantle  shear  veloci¬ 
ties  lower  than  previously  reported  values  for  this 
region. 


INTRODUCTION 


The  China  mainland  consists  of  a  number  of  tectonic 
subplates,  each  of  geologically  and  geophysically  distinct 
structures.  Based  on  an  investigation  of  the  surface 
geology,  Sun  and  Teng  (1977)  divide  China  into  4  basic 
tectonic  units  (Figure  1) :  (1)  the  Ching-Tibet  subplate, 

(2)  the  North-China  subplate,  (3)  the  Northwest-China 
subplate,  and  (4)  the  South-China  subplate.  The  present 
study  focuses  its  attention  on  the  South-China  subplate. 
Specifically  it  makes  use  of  a  group  of  closely  spaced 
earthquakes  which  give  a  number  of  pure-path  surface  wave 
trains  recorded  at  a  single  SRO  station.  The  resulting 
dispersion  curves  give  a  measure  of  the  observational  un¬ 
certainty  of  the  surface  wave  data. 

Using  a  non-linear  least  squares  procedure  based  on 
the  stochastic  model  (Franklin,  1970,  Jordan  and  Franklin, 
1970,  and  Derr  et  ^. ,  1970),  intermediate-  and  long- 
period  group  velocity  data  were  inverted  to  obtain  a 
detailed  shear  velocity  model  for  the  crustal  and  upper 
mantle  structure  of  the  South-China  subplate  (SOCHSP) . 
Recently  a  number  of  large  earthquakes  and  aftershocks 
have  occurred  at  the  boundary  between  the  Ching-Tibet 
subplate  and  the  South-China  subplate.  Those  events  have 
been  recorded  at  an  SRO  station  at  Taipei,  Taiwan  with 
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FIGURE  1.  TECTONIC  SUBPLATES  OF  CHINA 


the  -station  code  TATO.  SOCH£P  is  ideally  suited  for  a 
surface  wave  study  for  three  reasons:  (1)  the  South- 
China  subplate  can  be  considered  nearly  laterally  homo¬ 
geneous  llui-g,  1974;  Sun  and  Teng,  1977);  epicentral 
distances  are  adequate  to  allow  an  intermediate  to  long 
period  dispersive  wave  train  to  develop;  (3)  seismic 
events  occur  and  are  recorded  within  the  subplate,  thus 
eliminating  the  effects  of  crossing  plate  boundaries. 

Lacking  good  geophysical  data,  most  of  the  previous 
seismological  studies  were  either  limited  to  a  simple 
estimate  on  the  crustal  thicknesses  and  shallow  crustal 
structures  or  lacked  a  discussion  of  model  resolution  and 
standard  deviation.  Knowledge  about  the  crustal  e-nd  mantle 
shear  velocity  structure  of  SOCHSP  is  poor.  Tseng  and  Sung 
(1963)  determined  crustal  thickness  by  applying  the  tri¬ 
partite  method  (Press,  1956)  to  two  events  from  the  New 
Britain  Islands  for  phase  velocity  information.  The 
inversion  results  of  their  study  are  seen  in  Figure  2. 

The  numbers  inside  each  triangle  represent  the  average 
crustal  thickness  within  three  stations.  Using  short 
period  (4-12  seconds)  Rayleigh  and  Love  wave  group 
velocity  dispersion.  Sung  et  a^.  (1965)  determined  the 
thickness  of  the  sedimentary  layers  for  several  regions 
in  China.  Tung  (1974)  made  use  of  long  period  group 
velocity  data  and  tried  to  obtain  a  shear  velocity 
structure  for  SOCHSP.  But  the  method  he  used  was 
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essentially  a  trial-and-error  surface  wave  inversion 
technique  and  thus  he  was  unable  to  give  a  discussion 
about  resolution  or  standard  deviation  of  his  resulting 
models. 

Previous  work  has  illustrated  that  Rayleigh  waves 
produce  more  information  about  earth  models  than  Love 
waves  (Wiggin,  1972;  Derr  and  Landisman,  1972).  Wiggins 
(1972)  has  shown  that  by  simultaneously  inverting  Love 
and  Rayleigh  wave  data,  it  is  possible  to  obtain  both 
shear  wave  and  density  structures.  Though  group  velocity 
and  phase  velocity  contribute  almost  the  same  information 
(Bloch,  1969,  Bloch  et  a^. ,  1969,  Wiggins,  1972),  Bloch 
(1969)  and  Fix  (1975)  come  to  the  conclusion  that  group 
velocity  is  more  sensitive  to  the  earth  model  and  there¬ 
fore  should  produce  better  resolution  in  the  inversion. 
Through  a  multiple  filtering  technique,  this  study  first 
obtains  group  velocity  measurements  from  the  fundamental 
mode  Love  and  Rayleigh  waves,  followed  by  a  Gilbert- 
Backus  (1967,  1968)  type  inversion  to  obtain  a  crustal 
and  upper  mantle  structure  of  SOCHSP  together  with  an 
estimate  of  model  standard  deviation  and  resolution. 
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SOUTH-CHINA  SUBPLATE 

The  South-China  subplate  (SCS)  is  bounded  on  the 
south  and  east  by  the  Pacific  Ocean,  on  the  west  by 
longitude  103®E,  and  on  the  north  by  latitude  32 ®N 
(Figure  3) .  From  previous  gravity  and  seismic  studies 
(Tseng  and  Sun,  1963;  Tseng,  1973;  Tung,  1974)  the  crustal 
thickness  of  SCS  is  believed  to  increase  from  30  km  at 
the  eastern  boundary  to  50  km  at  the  western  boundary  with 
an  average  thickness  of  40  km. 

The  subplate  is  composed  of  two  parts,  the  Yangtze 
fault  block  in  the  northwest  and  the  South-China  folding 
block  in  the  southeast.  The  Yangtze  fault  block  is  a 
section  of  old  continental  crust  formed  during  the  Jinning 
Orogeny  (1300-950  m.y.),  and  was  subsequently  faulted  and 
strongly  metamorphosed.  Formed  during  the  Caledonian 
Orogeny  (570-400  m.y.),  the  South-China  folding  block 
was  a  section  of  transitional  (miogeosynclinal)  crust 
that  later  underwent  weak  metamorphism  and  folding. 
(Tectonic  Map  Compiling  Group,  1974,  Guizhou  Geologic 
Team  103,  1975) . 

While  surrounded  by  a  number  of  seismically  active 
zones,  SOCHSP  is  one  of  the  most  stable  regions  in  China. 
Though  a  number  of  large  events  occur  on  or  beyond  its 
western  boundary,  which  marks  the  India-Eurasia  collision. 


FIGURE  3.  THE  SOUTH-CHINA  SUDPLATE 


there  are  only  a  few  interplate  earthquakes  (Figure  4). 
According  to  Li  et  (1974) ,  this  can  be  explained  by 
the  fact  that  this  region  has  undergone  weaker  metamorphism 
than  other  regions  and  therefore,  the  rocks  are  softer  and 
do  not  lend  themselves  to  high  stress  accumulation. 

DATA 

The  data  used  in  this  study  are  surface  waves  from 
four  earthquakes  (Table  1)  focused  in  the  Szechwan 
Province  of  mainland  China.  The  events  were  recorded  at 
SRO  station  TATO.  A  number  of  events  over  one  path  were 
desired  in  order  to  establish  data  repeatability  and  to 
obtain  observational  uncertainties.  Because  of  instrument 
problems  and  poor  dispersion  results  (due  to  the  arrival 
of  wave  energy  at  certain  periods  that  have  taken  non- 
least-time  paths],  certain  components  of  ; ncividual 
earthquakes  were  not  used  in  obtaining  group  velocity 
information.  The  epicentral  region  and  path  are  shown 
in  Figure  3  (4-4 ' ) . 

SRO  data  in  digital  form  on  tape  is  convenient  for 
computer  processing  techniques  and  the  storage  of  a  large 
amount  of  information  on  a  single  data  tape.  The  SRO 
station  operates  at  high  sensitivity  thus  expanding  the 
detection  threshold  for  seismic  events  making  small 
magnitude  events  (M'v<5)  useful  for  surface  wave  analyses. 
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FIGURE 


TABLE  1 

LIST  OF  SEISMIC  EVENTS 


Date 

Focal 

Depth 

Epi central 
Distance 

Location 

(km) 

Magnitude 

(km) 

August  19,  1976 

32.893N 

104.189E 

33 

5.4 

1902.0 

August  22,  1976 

32.992N 

104.181E 

33 

4.7 

1884.0 

September  1,  1976 

32.460N 

104.152E 

18 

5.1 

1886.0 

September  3,  1976 

28.040N 

100.345E 

33 

5.2 

1854.0 

10 


With  this  new  system,  the  amount  of  data  available  for 
long-period  surface  waves  studies  is  greatly  increased. 

SRO  datatapes  contain  data  from  three  long-period 
channels,  sampled  at  the  rate  of  once  per  second  and  one 
vertical  component  sampled  at  the  rate  of  twenty  times 
per  second.  Long-period  data  is  recorded  continuously 
on  tape;  short-period  data  is  recorded  only  when  an  event 
has  been  detected,  and  is  Scunpled  at  the  rate  of  twenty 
times  a  second.  The  vertical,  north-south,  and  east-west 
long  period  channels  are  multiplexed  into  each  record. 

Each  SRO  data  record  consists  of  1000  words  (2000  bytes) . 
The  first  ten  words  are  header  information.  The  next  990 
words  are  data;  therefore,  there  are  5  minutes,  30  seconds 
of  data  per  long  period  record  (Peterson  et  al.,  1976). 

The  amplitude  and  phase  response  for  the  long-period 
instrument  is  given  in  Table  2.  An  analytic  form  of  the 
transfer  function  used  in  our  computer  program  is  found 
by  numerically  fitting  the  amplitude  and  phase  response 
into  a  polynomial.  The  transfer  function  has  the  form; 


n 

£ 

i 


ai  =  3.1635  X  lO"* 
ai  = 


.21163 


.14127 


wRumqMpmii 


j 

j 

1 

I 

I 

'1 

:3 


Frequency 

Relative  Amplitude 

Phase  Angle 
(Radians) 

.0667 

.460 

-1.33 

.05556 

.672 

-1.00 

.04546 

.912 

-0.63 

.04 

1.00 

-0.42 

.0333 

1.01 

-0.13 

.02778 

.9 

0.10 

.02326 

.724 

0.30 

• 

o 

N) 

.550 

0.47 

.01667 

.400 

0.63 

.0125 

.200 

0.92 

.01 

.108 

1.10 

.00667 

.031 

1.34 

.005 

.0104 

1.51 

.0025 

.00084 

1.84 

.00167 

.000184 

2.00 

. i 

_ .  . — - — 


bo  =  3.531  X  10-  bi  =  3.3316  x  10"^ 

b2  =  .12378  bs  =  1.0 

A  plot  of  this  transfer  function  is  given  in  Figure  5. 

Each  seismic  event  must  be  decoded  and  plotted  from 
the  data  tapes.  Pure  Love  waves  are  obtained  by  rotating 
the  data  coordinates.  An  example  of  a  plotted  seismic 
event  before  and  after  coordinate  rotation  is  given  in 
Figure  6  and  7. 

A  multiple  filtering  analysis  (Dziewonski  and  Hales, 
1970;  Hermann,  1973;  Tung,  1974;  Seekins  and  Teng,  1976) 
has  been  applied  to  the  surface  wave  train  to  obtain 
group  velocity  dispersion  curves.  A  typical  contoured 
result  of  this  type  of  analysis  is  shown  in  Figure  8. 
Because  of  the  epicentral  distances  involved,  only  periods 
of  up  to  65  seconds  were  used  for  Rayleigh  waves  and  5C 
seconds  for  Love  waves.  The  dispersion  curves  for  the 
fundamental  mode  Love  and  Rayleigh  waves  are  plotted 
in  Figures  9  and  10.  In  these  figures,  the  solid  line 
represents  the  average  group  velocities  obtained  from 
the  multiple  filtering  analysis.  The  statistical 
properties  of  the  observations  are  listed  in  Table  3  and  4. 


THEORY 

(loneralized  least  squares  methods  for  determining  the 
Eartli's  structure  from  seismic  data  have  been  formulated  by 
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3.44 
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TABLE  4 


LOVE  WAVE 
AVERAGED  OBSERVED 

PERIOD  GROUP  VELOCITY  (KM/SEC)  STANDARD  DEVIATION 


51.2 

3.95 

0.1 

42.667 

3.61 

0.1 

36.571 

3.52 

.03 

32.00 

3.41 

.03 

• 

28.444 

3.35 

.03 

25.60 

3.29 

.03  , 

23.273 

3.26 

.03 

21.333 

3.24 

.03 

19.692 

3.22 

.03 

18.286 

3.24 

.03 

17.067 

3.23 

.03 

16.00 

3.17 

.04 

15.059 

3.17 

.05 

> 

13.474 

3.17 

.03 

- 

12.19 

3.19 

.04 

11.13 

3.08 

.06 

!  • 

10.24 

3.09 

.03 

) 
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numerous  authors  (Jackson,  1972;  Wiggins,  1972;  and 
Crosson,  1976).  The  basic  logic  in  the  inversion  theory 
is  as  follows: 

The  observations,  Y^,  are  related  to  certain  model 
parameters,  X^,  in  some  known  way,  Y^  =  A(Xi ,X2 , . . • Xp) . 

A  quasi-linear  relationship  is  established  by  taking  a 
Taylor  series  expansion  about  some  ir^tial  value  Xj. 

.  3A.  , 

Y.  =  A. (X.)  +  — i  „iAX  +  higher  order  terms 

^  ^  *  3Xj  5Ce 

Yi  =  A. (X;)  +  AY. 

Ignoring  the  higher  order  terms,  we  obtain  the 
functional  relationship; 


AY^  = 


3Ai 

3Xj 


(1) 


where  AY^  =  observed  data  minus  calculated  observations. 
Equation  1  in  matrix  form  becomes; 


AAX  =  AY 


Where  A  is  a  matrix  whose  elements  are  the  partial 
derivatives  in  Equation  1,  AX  is  a  p  x  1  vector  whose 
elements  are  correction  to  the  initial  model  vector  xi, 
and  AY  is  the  difference  between  the  observed  and 
theoretical  data. 


In  our  case  the  observations  are  group  velocities  of 


! 

i 

1 


Love  and  Rayleigh  waves,  and  the  model  parameters  are 
shear  velocities  and  densities  for  a  horizontally  layered 
earth.  The  partial  derivatives  in  equation  1  are  generated 
numerically  by  the  method  described  by  Rodi  et  (1975) . 

The  classic  least  squares  solution  to  equation  2  is 
by  minimizing  the  euclidean  length  of  AAX-AY  or 

II  AAX-AY  II  -  (AAX-AY)'^(AAX-AY)  =  e 

Given  equation  2 ,  the  solution  is  (Hanson  and 
Lawson,  1975) : 

AAX  =  AY 

T  T 

A^AAX  *  A  AY 

^  T  1  T 

AX  »  (A'‘a)"'A^AY 

where  a"x  is  the  estimated  correction  vector. 

A  suitable  inverse  to  this  problem  can  be  obtained 
by  using  the  sing”lar  value  decomposition  (SVD)  of  the 
matrix  A. 

A  =  UAv"^ 

where  the  columns  of  U  are  the  eigenvectors  associated 
with  the  colums  of  A,  the  rows  of  V  are  the  eigenvectors 
associated  with  the  rows  of  A,  and  A  is  a  diagonal  matrix 
of  non-zero  eigenvalues  of  the  matrix  A. 
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T 

A  «  A^A  =  , 

• 

I 

1  0 

If  we  let 

h  =  (A  A)  A 

then 

HAAX  =  HAY 

and 

AX  =  HA  X 

and 

£x  =  HAY 

Lsing  the  SVD  we  obtain  (Jackson,  1972) : 

a"x  =  (AV'u'^)AY 
H  =  AV'U*^ 


(3a) 

(3b) 

(3c) 


{3d) 


Here,  H  is  known  as  the  pseudo-inverse  of  equation  2. 

The  matrix  product  HA=R  is  known  as  the  resolution 
matrix,  and  is  a  measure  of  the  uniqueness  of  the  solution 
(Jackson,  1972) .  The  rows  of  R  are  called  the  resolving 

kernels.  From  equation  3c,  it  can  be  seen  that  an  element 

^  ^  tih 

of  ax,  AXj^,  may  be  interpreted  as  convolving  the  v  row 

of  R  with  the  vector  AX.  Therefore,  AX^  can  be  considered 

the  weighted  sum  of  nearby  values  (Jackson,  1972) . 
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T 

Because  A  A  is  nearly  singular,  a  number  of  problems 
arise:  (1)  The  solution  vector  becomes  large  which  can 

cause  the  problem  to  leave  the  region  of  linearity;  (2) 
the  solution  oscillates  with  each  iteration.  By  looking 
at  the  pseudo-inverse  H  we  can  see  the  cause  of  this 
instability. 

H  =  VA-‘U^ 

where  the  matrix  A“‘  equals 


T 

and  the  X^'s  are  the  eigenvalues  of  the  matrix  A  A.  If 
T 

A  A  is  near  singular,  one  or  more  of  the  eigenvalues  will 
be  approaching  zero.  From  equation  3  it  can  be  seen  that 
a  small  eigenvalue  will  cause  a  large  change  in  one  or  more 

ys 

values  of  the  correction  vector  AX.  Similarily,  it  has 
been  shown  that  the  variance  of  the  model  parameters  are 
inversely  proportional  to  X  (Jackson,  1972) .  Thus  small 
X's  produce  large  standard  deviations. 

There  are  two  main  approaches  to  stabilizing  the 


inversion  process.  One  method  is  to  examine  the 
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eigenvalue  spectrum  of  the  matrix  A.  Small  eigenvalues 
are  removed  when  the  variance  becomes  too  large.  But 
removing  the  eigenvalues  degrades  the  resolution.  This 
trade-off  between  resolution  and  variance  has  been 
discussed  by  a  number  of  authors  (Jackson,  1972;  Wiggins, 
1972;  and  Crosson,  1976) .  Braille  and  Keller  (1975)  used 
this  method  for  the  inversion  of  group  velocity  data. 

Inversion  stability  can  also  be  achieved  by  using 
Marquardt's  method  (Marquardt,  1963)  or  the  stochastic 
inverse  (Franklin,  1970) .  These  procedures  lead  to  the 
suppression  of  small  eigenvalues.  Franklin's  stochastic 
inverse  estimates  the  parameter  in  the  presence  of  noise. 
Equation  2  becomes: 

AAX  +  n  =  AY  (4) 

where  n  is  a  vector  of  observational  noise.  Franklin 
(1570)  shows  that  the  solution  to  equation  4  is  given  by: 

AX  =  Wa'^(AWa'^  +  ED)"‘  AY  (5) 

where  W  is  the  covariance  matrix  of  the  parameter  and  ED 
is  the  covariance  matrix  of  the  observations. 

Marquardt's  method  considers  a  minimization  of  the 
functional 

(AAX  -  AY)'^(Ed'^ED)  (AAX  -  AY)  -  Ax'^(w'^W)  AX  (6) 
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with  the  solution  given  by: 


AX  =  a  +  Wp-')“*A'^ED'*AY  (7) 

T  I  T 

Letting  A  LD"*  =  A  ,  equation  7  becomes: 

AX  =  (a‘^A  +  0^1)  a'^AY  (3) 

where  o*  is  the  variance  of  the  parameter. 

By  looking  at  the  SVD  of  equation  8  we  can  see  how 
Marquardt's  method  tapers  the  eigenvalue  spectrum; 

AX  =  V|  (A^  +  o^D'^Alu'^AY  (9) 

An  element  in  brackets  from  equation  9  equals;  (Crosson, 
1976) 


I  (A^  +  a*I)*‘A|  .  =  X./(X2^  +  oM 

So  as  X^  goes  to  0,  the  parameter  AX^  goes  to  zero.  A 
similar  analysis  for  the  stochastic  inverse  shows  that  it 
too  suppresses  small  eigenvalues. 

The  quantity  can  also  be  thought  of  as  a  trade¬ 
off  parameter  between  resolution  and  variance  (Der  ^ 
al. ,  1970;  Wiggins,  1972;  Crosson,  1976).  In  Marquardt's 
method,  the  resolution  equals  (Crosson,  1976}  : 

R  »  HA  »  (a'^A  +  o*I)A*^A  *  V|  (A  +  0*1)"' A^  IV*^  (10) 


For  the  stochastic  inverse 


(11) 


R  »  I  (A^  +  jW*^ 

The  quantity  in  brackets  of  equations  10  and  11  is  a 
diagonal  matrix  whose  elements  have  the  form: 

If  0*  =  0,  then  R  =  I,  but  as  o*  goes  to  0  the  variance 
of  the  par^uneters  becomes  larger.  So  we  adjust  the  size 
of  0*  until  there  is  an  acceptable  trade-off  between 
standard  deviation  and  resolution. 

Two  other  elements  that  must  be  analyzed  in  the 
inversion  process  are  the  apriori  (SDXo)  and  postpriori 
(SOX)  standard  deviations  of  the  inversion  parameters; 

SDXo  =  d'^(oWP"‘)D 

SDX  =  /d'^(a'^ED"‘ A  +  oWP”M"'D 

where  D  =  a  row  vector  of  a  delta  matrix  =  (0,0,...,1,...0). 
This  analysis  will  give  a  measure  of  how  much  new  informa¬ 
tion  is  coming  from  the  data.  If  SDXp  and  SDX  are 
approximately  the  seune,  we  know  that  very  little  informa¬ 
tion  is  coming  from  the  data. 

The  advantages  of  Marquardt's  method  and  the  stochas¬ 
tic  inverse  is  that  the  SVD  of  A  is  not  explicitly  deter¬ 
mined  and  there  is  no  decision  to  be  made  on  the  rank  of 
T 

the  matrix  A  A. 

A  combination  of  these  two  methods  are  used  in  our 
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inversion  process.  When  the  number  of  parameters  is 
greater  than  the  number  of  data,  the  stochastic  inverse 
is  used.  If  the  number  of  observations  is  greater  than 
the  number  of  parameters,  Marquardt's  method  is  used. 
This  procedure  leads  to  the  minimum  number  of 
computations. 


DISCUSSION 

Initial  models  for  the  inversion  were  based  on  the 
results  of  a  previous  surface  wave  study  done  by  Tung 
and  Teng  (1974).  Their  final  model  was  obtained  from 
Love  and  Rayleigh  wave  group  velocity  data.  Our 
starting  model  (Ml,  Table  5)  was  inverted  for  both  shear 
velocity  and  density. 

During  the  inversion,  large  instabilities  of  the 
density  parameter  were  found  in  the  crust.  This  was 
corrected  by  changing  the  starting  model  and  inverting 
for  compressional  wave  velocity.  In  the  derived  model 
(DMl)  a  low  velocity  layer  (LVL)  can  be  inferred  in  the 
sixth  layer,  but  the  validity  of  this  LVL  is  questionable 
because  of  the  insufficient  data  in  our  study.  Examin¬ 
ation  of  the  apriori  and  postpriori  standard  deviations 
for  the  shear  velocity  parameter  (Table  5)  shows  that 
there  is  only  a  small  amount  of  information  in  the  data 
applicable  to  depths  below  the  fifth  layer.  Thus  longer 
period  information  is  needed  to  conJirm  velocities  beyond 
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table  5 


STARTING  MODEL 
Ml 


LAYER 


THICKNESS 

(KM) 

D-VELOCITY 

(KM/SEC) 

SHEAR  VELOCITY 
(KM/SEC) 

DENSITY 

10 

5.0 

3.0 

3.0 

10 

5.0 

3.1 

3.05 

20 

6.0 

3.5 

3.2 

20 

8.40 

4.6 

3.3 

25 

8.05 

4.55 

3.34 

50 

7.80 

4.2 

3.37 

50 

8.00 

4.3 

3.37 

FINAL  MODEL 
DM1 

SHEAR  VELOCITY 

(KM/SEC)  SDXo  SDX 


3.245 

.277 

.046 

3.193 

.277 

.076 

3.767 

.196 

.055 

4.498 

.196 

.106 

4.495 

.175 

.112 

4.283 

.124 

.105 

4.435 


124 


120 


this  depth. 

A  new  initial  model  (M2)  using  thinner  layers  in 
the  crust  and  mantle  was  inverted  for  shear  wave  velocity, 
density  and  compressional  wave  velocity.  The  mantle 
shear  velocities  for  M2  were  also  modified,  as  were 
density  and  P  velocity  values  in  the  crust  (Table  6) . 

By  examining  the  relative  sizes  of  the  partial 
derivatives  of  different  parameters,  we  can  determine 
the  importance  of  each  parameter  in  the  inversion.  The 
partial  derivatives  of  group  velocity  with  respect  to 
P  wave  velocity,  shear  velocity  and  density  for  the 
fundamental  Rayleigh  wave  mode  are  plotted  in  Figure  11a, 
b,  and  c.  The  partial  derivatives  of  group  velocity  with 
respect  to  shear  velocity  and  density  are  plotted  in 
®’igure  12a,  b  for  the  fundamental  Love  wave  mode.  Figures 
11  and  12  show  a  number  of  important  factors  that  affect 
this  study: 

(1)  The  partial  derivatives  for  P  wave  velocity  in 
the  upper  crust  are  not  insignificant  when  compared 
to  density  ^md  shear  velocity  and  therefore,  affect 
the  inversion. 

(2)  The  P  wave  partial  derivatives  in  the  lower 
crust  and  mantle  arc  small,  and  are  not  significant 
in  the  inversion. 

(3)  The  partial  derivatives  with  respect  to  shear 
velocity  are  larger  than  those  for  density  and  P 
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TABLE  6 


STARTING  MODEL 
M2 


LAYER 

THICKNESS 

P-VELOCITY 

(KM/SEC) 

SHEAR  VELOCITY 
(KM/SEC) 

DENSITY 

5 

5 . 5 

3.0 

2.6 

5 

5.5 

3.0 

2.8 

5 

3 . 5 

3.1 

3.05 

5 

5.5 

3.1 

3.05 

10 

6.0 

3.5 

3.1 

10 

6.0 

3.5 

3.1 

10 

8.1 

4.6 

3.3 

10 

8.1 

4.6 

3.3 

10 

8.0 

4.4 

3.3 

lb 

8.0 

4.4 

3.3 

25 

3.0 

4.4 

3.37 

25 

8.0 

4.4 

3.37 
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FIGURE  12a 


wave  velocity.  Therefore,  the  resolution  for  shear 
velocities  is  better  than  for  the  two  other 
parameters. 

The  starting  model  M2  and  the  derived  model  (SCSI) 
with  standard  deviation  error  bars  are  shown  in  Figure  13 . 
The  fit  of  theoretical  group  velocities  for  SCSI  are 
shown  in  Figures  14  and  15.  The  averaging  Jcernals  are 
listed  in  Table  7,  where  only  the  elements  of  the  kernels 
corresponding  to  the  shear  parameters  are  included.  The 
other  elements  for  P  wave  velocity  and  density  are  small 
and  do  not  affect  the  structure.  The  width  of  these 
kernels  indicate  that  it  is  possible  to  resolve  layers 
to  a  thickness  of  5  km  in  the  upper  crust,  10  km  in  the 
middle  crust  and  20  km  in  the  lower  crust.  However, 
resolution  in  the  mantle  is  quite  poor.  All  the  standard 
deviations  were  found  to  be  less  than  .08  km/sec. 

Increasing  the  resolution  in  the  lower  crust  and  mantle 
(by  using  a  different  trade-off  parameter  a)  leads  to 
instability  in  the  shear  velocity  model,  and  unacceptably 
large  standard  deviations.  These  features  could  be 
improved  if  higher  mode  data  were  available. 

The  results  of  the  inversion  indicate  that  only  a 
few  layers  are  justified  to  describe  the  crust  with  the 
given  dispersion  information.  As  an  option  in  the  inver¬ 
sion  program,  averaging  kernals  can  be  made  to  resemble 
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a  desired  box  car  vector.  An  exaunple  would  be  if  we 
wished  to  average  the  first  two  parameters,  then  we  would 
want  the  kernel  to  resemble  a  box  car  function  such  as 
described  in  Table  8a.  In  the  inversion,  shear  wave 
parameters  from  layers  2,  3  and  4  were  averaged,  as  were 
parameters  from  layers  3  and  6,  and  layers  7  and  8, 

Table  8c.  The  actual  kernels  for  these  averaged  para¬ 
meters  are  shown  in  Table  8b.  The  averaged  values  and 
standard  deviations  are  plotted  in  Figure  16.  (Model 
ASCSl,  Table  9).  A  comparison  with  Tung  and  Teng  (1974) 
results  indicates  here  that  crustal  shear  velocities 
are  higher  than  previously  reported  values.  And  the  lid 
velocity  in  the  mantle  is  found  to  be  slightly  less  than 
previously  determined.  A  LVL  at  GO  km  is  consistent  with 
our  dispersion  data. 

A  list  of  the  apriori  and  postpriori  standard 
deviations  for  SCSI  is  given  in  Table  10.  There  is  almost 
no  new  information  that  can  be  derived  from  the  data  past 
70  km.  In  the  mantle  in  general,  the  amount  of  new 
information  from  the  data  is  quite  small.  But  the 
persistence  of  the  LVL  at  the  60-70  km  depth  might  indicate 
a  true  velocity  reversal  and  possibly  the  transition 
between  the  lithosphere  and  asthenosphere . 
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FIGURE  16 


TABLE  9 


LAYER 

iICKNESS 

(KM) 

MODEL  ASCS  1 

SHEAR 

VELOCITY 

(KM/SEC) 

SDXo 

SDX 

' 

5 

3.03 

.347 

.087 

15 

3.36 

.200 

.062 

20 

3.84 

.173 

.063 

20 

4.51 

.173 

.120 

10 

4.28 

.245 

.210 

15 

4.36 

.200 

.175 

25 

4.44 

.155 

.142 

25 

4.44 

.155 

.150 
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TABLE 

FINAL  MODEL 
SCPl 


SHEAR  VELOCITY 

SDXo 

SDX 

3.03 

.347 

.087 

3.43 

.347 

.169 

3.24 

.347 

.231 

3.4 

.347 

.234 

3.89 

.245 

.150 

3.78 

.245 

.169 

4.49 

.245 

.203 

4.53 

.245 

.210 

4.28 

.200 

.175 

4.36 

.155 

.142 

4.35 

.155 

.150 

47 


CONCLUSIONS 


By  comparing  the  size  of  the  P-wave  partial 
derivatives  with  density  and  shear  wave  partial  deriva¬ 
tives,  it  was  determined  that  compressional  wave  velocity 
cannot  be  neglected  when  inverting  surface  wave  data. 

The  resolution  for  shear  velocity  is  better  than  the 
resolution  for  density  and  p-wave  velocity. 

The  inversion  of  multi-layer  models  indicates  that 
a  3  layer,  40  km  crust  adequately  describes  the  South- 
China  subplate.  The  resolution  ranged  from  5  to  20  km  in 
the  crust  with  model  standard  deviations  of  less  than 
.08  km/sec.  The  shear  velocities  were  found  to  be  higher 
than  those  velocities  reported  by  Tung  (1974) . 

Upper  mantle  shear  velocities  were  found  to  be  lower 
than  previously  reported  values.  Resolution  in  the  mantle 
was  quite  poor  as  no  dispersion  information  is  obtained 
from  this  study  beyond  65  seconds.  The  apriori  and 
postpriori  standard  deviation  indicated  that  very  little 
new  information  was  being  obtained  from  the  data  about 
the  mantle.  But  because  of  the  persistence  of  a  LVL 
between  60  and  70  km,  it  is  presumed  that  transition 
between  the  lithosphere  and  asthenosphere  occurs  here. 
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APPENDICES 


In  tis  appendix  attention  will  be  focused  primarily 
on  how  to  use  the  computer  programs  discussed  in  this 
thesis.  There  are  very  little  mathematical  and  algorith¬ 
mic  descriptions  given  in  this  text.  These  descriptions 
are  covered  in  the  reference  section  in  Table  A.l,  and  it 
is  to  the  benefit  of  the  user  to  read  these  papers  before 
attempting  to  use  the  programs.  The  algorithms  have 
been  tested  and  implemented  in  Fortran  4,  with  the  entire 
set  of  codes  presented  in  card  image  format  or  set  into 
a  load  module. 

The  3  main  programs,  PLOT,  GVA  (group  velocity 
analysis)  and  SWI  (surface  wave  inversion) ,  are  presented 
here  in  the  order  used  in  the  actual  surface  wave 
inversion  process. 


PLOT 

After  obtaining  the  SRO  data  tapes,  the  problem  is 
to  put  the  data  into  a  recognizable  form;  this  is 
accomplished  by  the  program  PLOT. 

A.  Deck  Setup. 

Figure  A.l  shows  the  order  of  cards  in  the  deck. 
The  only  cards  that  need  to  be  changed  are  the  input 
cards  and  the  "VOL*SER“"  card.  The  "VOL*SER=kXX" 
number  is  simply  the  storage  location  or  "bin" 
number  at  the  University  Computer  Center.  It  tells 
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TABLE  A.l 


PROGRAM 

PLOT 

GVA 

SWI 


DESCRIPTION  REFERENCE 

Rotates  and  plots  - 

seismogram 

Group  velocity  Tung  (1974)  and 

analysis  of  rotated  Hermann  (1972) 

seismogr^un 

Does  actual  surface  Wiggins  (1972) 

wave  inversion  to  and  Jackson  (1972) 

obtain  theoretical 
earth  model 


INPUT  DATA  DECK 


the  operator  where  to  find  your  tape. 

B.  Input  Cards. 

Input  cards  and  their  formats  are  listed  below. 

Card  1  FORMAT  (F5.2) 

THETA  “  the  angle  coordinate  axes 
rotation 

Card  2  FORMAT  (2(I4,6X)) 

IFIRST  =  the  value  of  the  first  SRO  data 
record  to  be  read  in 
NDATA  =  number  of  data  points  to  be 
plotted  (maximum  =  2000,  if 
ISKIP  =  0) 

Card  3  FORMAT  (2(I4,6X) 

ISKIP  =  the  number  of  short-period 

records  that  are  skipped  to  get 
to  the  next  long-period  record 
NEXTRA  =  the  number  of  extra  data  points 

to  be  plotted.  There  can  be  more 
than,  but  NEXTRA  +  NDATA  cannot 
be  greater  than  2000. 

Card  4  FORMAT  (2(I4,6X)) 

ISTRT  =  the  point  at  which  the  user  wants 
punched  output  to  begin 
PNCH  »  the  number  of  points  to  be  punched 
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C.  Output. 

The  output  consists  of  three  plots  (for  an 
example  see  Figure  3) .  If  punched  data  is 
requested  then  output  will  have  the  format 
(5 (I4 ,E12 .6) ) .  The  vertical,  north-south,  and 
east-west  component  punched  output  are 
separated  by  a  bl2m)c  card  and  a  card  that  says 
"NUM  =  . "  The  punched  output  is  ready  to  use 

in  the  multiple-filtering  program,  GVA. 


GVA 

The  multiple-filtering  analysis  of  the  digitized 
surface  wave  train,  is  accomplished  in  the  program  GVA. 

A.  Deck  Setup. 

Figure  A. 2  shows  the  order  of  cards  in  the  deck. 
The  only  cards  that  need  to  be  changed  are  the 
input  cards. 

B.  Input  Cards. 

Input  cards  and  their  format  are  listed  below. 
Card  1  FOWiAT  (20A4) 

NAME  =  an  alphanumeric  description  of 
the  data 

Card  2  FORMAT  (15,  fS.O) 

LX  =  the  total  number  of  data  points 

to  read  in  (this  value  has  to 
equal  2**N 
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DX  =  The  digitization  interval  (for 

SRO  long-period  data  this  value 
equals  1) 

Card  3  FORMAT  (2F10.3) 

DIST  =  the  epicentral  distance 

TCON  =  the  difference  between  the  earth¬ 
quake  origin  time  and  the  time  at 
the  start  of  the  first  record 
using  in  program  PLOT 

Card  4 

X(I)  =  punch  output  of  PLOT 
Y(I)  =  punch  output  of  PLOT 

C.  Output 

The  most  important  part  of  the  output  is  a  plot 
of  the  group  velocity  vs.  period.  (For  an 
example  see  Figure  5.) 

D.  Problems  using  GVA 

1.  Make  sure  there  are  2**n  data  points  (where 
n  is  some  integer  value) 

2.  Have  punched  output  from  that  part  of  the 
seismogram  which  represents  the  group 
velocity  arrivals  between  4.5  km/sec  and 
2.0  km/sec. 

3.  You  might  have  higher  mode  arrivals  in  your 
dispersion  curve.  This  is  not  a  problem. 

Try  to  identify  these  arrivals  and  use  them 
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INPUT  DATA 


4. 


in  your  inversion. 

If  you  have  plotted  the  seismogram,  (p  wave, 
s  wave,  and  surface  wave  arrivals)  TCON 
might  be  negative.  Input  this  the  same  way 
you  would  for  any  F  type  format. 

SWI 

SWI  takes  the  group  velocity  data  obtained  in  GVA  and 
inverts  it  for  an  earth  model.  The  program  has  been  set 
up  into  a  load  module.  Therefore,  the  only  part  of  the 
program  that  has  to  be  submitted  to  the  computer  is  the 
input  data  deck  and  some  job  control  language  cards. 

Figure  A. 3. 

The  comment  cards  in  the  listing  of  SWI  give  an 
excellent  description  of  the  different  input  and  output 
parameters  (Appendix  3).  This  listirg  also  shows  how  to 
input  the  data  into  the  load  module.  The  discussion  here 
will  be  limited  to  the  problems  in  using  the  program. 

A.  Problems  in  using  SWI 

1.  Invert  for  shear  velocity,  and  density.  If 
you  have  a  problem  in  obtaining  a  reasonable 
final  model,  then  invert  for  compressional 
velocity  also. 

2.  When  choosing  values  for  "TOP"  use  values 
like  10“*,  10"*,  10"*,  and  lO"**.  Do  not 
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use  values  in  between  these.  Probably  the 
smallest  value  you  will  need  to  try  is  10'** 
I  found  that  3  or  4  iterations  at  10“ *  got 
the  necessary  RMSFIT. 

You  want  your  final  RMSFIT  at  less  than  1. 
But  anything  up  to  a  RMSFIT  of  1.5  is 
probably  reasonable. 

Remember  we  are  trying  to  get  stable 
models.  St2d}ility  is  achieved  when  all 
your  PNMPG  are  less  than  STDEV  K. 


Ir  this  section  an  illustration  is  given  on  how  to 
use  the  programs. 


EXAMPLES 

An  earthquake  occurred  in  the  Szechwan  Province  of 
mainland  China  and  was  recorded  at  station  TATO.  The 
origin  time  was  9/1/76,  1:6:51.8,  and  the  epicentral 
distance  is  1886.28  km. 


PLOT 

The  angle  of  axis  rotation  was  found  to  be  21.7®.  Ws 
now  want  to  know  when  the  surface  wave  arrives  at  TATO. 
Divide  the  epicentral  distance  by  4.5  km/sec  and  add  this 
value  to  the  origin  time.  This  value  equals  1:13:27.6. 

Now  find  out  when  this  time’  occurs  in  the  SRC  datatape. 

Jr.  cur  example,  the  first  arrival  occurs  in  record  240. 

GVA 

We  now  have  punched  output  from  PLOT  and  wish  to  get 
a  dispersion  curve.  We  know  that  the  epicentral  distance 
is  1886.28.  To  get  TCON  subtract  the  origin  time  by  the 
starting  time  of  record  240;  this  value  equals  336.1. 

It  should  be  noted  that  I  only  have  punched  output  from 
150-450  seconds,  a  total  of  301  points.  But  I  need  512 
points  (2**9) .  This  problem  is  remedied  by  adding  on  the 


66 


end  of  the  punched  output  (from  451-661)  a  set  of  points 
whose  amplitude  values  equal  zero. 


r 


> 


1 
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In  this  section  the  input  data  sequence  is  given  for 
the  surface  wave  inversion  program  {Rodi,  1977)  . 

A.  Input  Cards 

Card  1  FORMAT  (812, 2E8 . 0. 3F8 . 0) 

LEAST  =  Flag  controlling  how  calculations 
are  done.  Normally  use  LEAST  =  0 
for  NrARK>NDATA,  LEAST  =  1 
otherwise 

NITERX(l)  Number  of  iterations  using  ALF 
TOP  on  this  card 

NALF  Number  of  extra  ALF/TOP  pairs 

read  in  later 

IFVX  FLAG  IFVX  *  1  Calculate  and  print 

out  'St  Dev  x' 

=  0  Don't  do  it 

IFVK  FLAG  IFVK  =  1  Calculate  and  print 

out  'St  Dev  k' 

*  0  Don't  do  it 

IFKERN  FLAG  IFKERN*  2  Print  and  punch 

averaging  kernels 

*  1  Print  averaging 

kernels 

a>  0  Don't  do  it. 

NJUMP  >  Number  of  numbers  on  card  12, 

which  describes  which  parameters 
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not  to  calculate  variances  for 


N0RI4  FLAG  WORM  »  1  Normalize  averaging  ^ 

kernals  read  in 
later 

«  0  Take  averaging 
kernals  read  in 
later  as  is 

TOPX(l)  Value  of  TOP.  If  TOP>0  it  is  used 
to  compute  ALF 

ALFX(l)  Value  of  ALF 

QUITl  ,  Stopper.  It  stops  iterating  when 
RMSFIT<QUIT1 

QCJIT2  Stopper.  It  stops  iterating  when 

:<MSPIT (present)  >Quit  2*RMSFIT  f 

(last) 

THRESH  Print/punch  only  the  averaging 

kernals  for  which  'QK'<THRESH 
If  NALF>0,  then  read  in  .NALF.  cards  of  Card  2 
Card  2  FORMAT  (2X, 12, 12X, 2E8. 0) 

NITERX  Number  of  iterations  using  ALF/ 

TOP  on  this  card.  If  NITERX  *  0, 
then  (1)  subsequent  Card  2's  must 
have  NITREX  “0,  and  (2)  parameter 
corrections  gotten  using  this. 

ALF/TOP  will  not  be  used  to  get 
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final  model.  NITERX=0  is  used  to 
get  RX-STDEVX  and  qk-ETDEVK 
trade-off  curves  for  next-to-last 
model. 

TOPX  Value  of  TOP 

ALFX  Value  of  ALF 

Card  3  FORMAT  (511,13) 

IFLATR  FLAG  IFLATR=1  Do  sphericity 

correction  for 
Rayleigh  using 
'Biswas  Earth 
flattening ' 

=0  Don't  do  it 

IFLATL  FLAG  IFLATL«1  Do  sphericity 

correction  for 
Love  using  "Biswas 
Earth  flattening' 
=0  Don't 

IFDIS  FLAG  IFDIS  =1  Print  out  the 

layer  displace¬ 
ments 

=0  Don't  do  it 

IFENG  FLAG  IFENG  »1  Print  out  the 

layer  energy 
■0  Don't  do  it 

IFPAR  FLAG  IFPAR  *1  Print  out  the 
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layer  partial 
derivatives  W.R.T. 


KEEP 


KEEP=0: 
Note ; 


Card  4 
ID 

Card  5 


NDATA 

LOGO 


layer  parameters 
*0  Don't  do  it 
Plus  or  minus  a  data  set  reference 
no.  {say  a  tape  unit)  for  saving 
inversion  results.  For  example: 
KEEP=91:  Save  inversion  results 
as  of  last  iteration 
KEEP--91:  Save  this  and  data  and 

partials  for  final  model 
Save  nothing 
Data  and  partials  for 
final  model  are  only 
computed  if  they  are 
requested  to  be  saved 

FORMAT  (18A4) 

Alphanumeric  data  descriptor 
FORMAT  (I4,2I1,2F8.0) 

Number  of  data. 

FLAG  LOGD=l  Inverse  program 

uses  natural  log 
of  C  &  U  for  obser¬ 
vations,  so  read  in 
standard  devs  of 
log  C  or  log  U 
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=0  C  and  U  are  used 
for  observations 

IFTORF  FLAG  IFTORFS2  Velocities  are 

functions  of  period 
=3  Velocities  are 
functions  of 
frequency 

DTFLOG  Pertubation  of  log  frequency  used 

to  compute  U-partials.  .001  is  a 
good  value 

STEP  Trail  phase  velocity  increment/ 

decrement 

Number  of  card  6's  =  number  of  frequencies 
Number  of  card  7's  »  number  of  data 
Card  6  FORMAT  (4I1,2F8.0) 

KINDX  FLAG  KINDX=1  For  Rayleigh  Wave 

-2  For  Love  Wave 

MODEX  Mode  number,  i.e.,  0  for  the 

fundamental  mode,  1  for  the 
first  higher  mode  and  so  on. 

IFCUX(l)  FLAG  IFCUX(1)«1  Read  in  phase 

velocity. 

-0  No. 

IFC0X{2)  FLAG  IFCUX{2)«1  Read  in  group 

velocity. 

«0  No. 


73 


IF  Period  for  IFT0RF=2  or  frequency 

for  IFTORF-3. 

Start  Guess  root  (Phase  velocity)  for 
the  period/frequency. 

There  follow  .IFCUX(l)  +  IFCUX(2). 

Cards  of  Card  7 

Card  7  FORMAT  {F8. 0, 14 , 5 (14 ,F8 . 0) / (6 (14 . FO. 0) ) 
DOBS  Observed  value  of  phase (C)  or 

group (U)  velocity 

NVAR  Nxunber  of  data  this  obs  has 

covariances  with  (including  itself 
but  not  data  already  read  in) 

MAP  Datum  numbe*  this  obs  covaries 

with  minus  datum  number  of  this 
obs  (thus  Map“0  means  variance# 

ED  Value  of  standard  deviation  (co- 

deviation)  .  If  MAP>0:  Codev  = 
Sign(covar) *SQURT(Abs (Covar) ) . 

Card  8  FORMAT  (18A4) 

ID  Alpnanumeric  model  descriptor. 

Card  9  FORMAT ( 14 , 411) 

Layers  Number  of  layers  including  water 
layer  and  half  space  in  the  guess 
model. 

LOGP  FLAG  LOGP  >1  Inverse  prog  uses 
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log  VP,  VS,  density 

as  parameters. 

=0 

VP,  VS,  Density 

used  as  parameters. 

IFVP 

FLAG 

IFVP 

=1 

Invert  for  P-wave 

velocity  (VP) 

=0 

No. 

IFVS 

FLAG 

IFVS 

=1 

Invert  for  S-wave 

velocity  (VS) 

«o 

No. 

IFDEN 

FLAG 

IFDEN 

*1 

Invert  for  density 

=0 

No. 

There  now  follow  .layers,  cards  of  card  10 
Card  10  FORMAT  (4F10.4) 

HQ  Layer  thickness  in  km. 

?.LPQ  Layer  P-V’ave  velocity  in  km/sec. 

BETQ  Layer  S-Wave  velocity  in  km/sec. 

RHOQ  Layer  density  in  gm/cc. 

There  now  follow  . IFVP  +  IFVS  +  IFDEN . 

Sets  of  card  11.  Number  of  cards  in  each 
set  depends  on  how  the  user  divides  the 
model  into  segments. 

Card  11  FORMAT  (I4,2F8.0) 

T.AY2  Number  of  bottom  layer  of  segment 

using  WIDTH  and  DEV  given  on  this 
card. 
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WIDTH 


Smoothing  interval  in  km. 
Recommendation:  Set  WIDTH  =  0. 

DEV  If  WIDTH  =  0,  DEV**2/ALF  is 

assumed  variance  per  km  of 
pareuneters  in  this  segment. 

Note:  Absolute  size  of  DEV  is 
unimportant.  Relative  sizes  of 
DEV  in  different  segments  allows 
constraining  of  some  parameters 
compared  to  others.  (Small 
DEV=>Big  constrain.)  Absolute 
degree  of  constraining  can  be 
controlled  with  ALF.  (Big  ALF  = 
>Big  Constrain.) 

If  NJUMP  not  equal  to  zero,  use  card  12 
Card  12  Free  FORIIAT  (Use  only  columns  1-72) 

On  this  card  are  NJUMP  Numbers  which  indi¬ 
cate  the  parameters  for  which  variances 
and  averaging  kernels  are  not  to  be  calcu¬ 
lated.  For  exeunple,  if  NJUMP  »  4  and  this 
card  contains  1,  -22,  25,  and  +30,  then 
they  won't  be  done  for  parameters  1  through 
22,  parameter  25  and  parameter  30.  The  '+' 
sign  before  30  can  be  omitted. 

Card  13  Free  FORMAT  (Use  only  columns  1-72) 

This  card  contains  an  ideal  averaging  kernal 
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(the  vector  X).  Example:  4*0  10*. 1 

999*0.  This  vector  averages  pareuneters 
through  14. 


